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Abstract 

In this work, a stochastic, physics-based model for Lithium-ion batteries (LIBs) 
is presented in order to study the effects of parametric model uncertainties on the 
cell capacity, voltage, and concentrations. To this end, the proposed uncertainty 
quantification (UQ) approach, based on sparse polynomial chaos expansions, re¬ 
lies on a small number of battery simulations. Within this UQ framework, the 
identification of most important uncertainty sources is achieved by performing 
a global sensitivity analysis via computing the so-called Sobol’ indices. Such 
information aids in designing more efficient and targeted quality control proce¬ 
dures, which consequently may result in reducing the LIB production cost. An 
LiC6/LiCo02 cell with 19 uncertain parameters discharged at 0.25C, 1C and 4C 
rates is considered to study the performance and accuracy of the proposed UQ 
approach. The results suggest that, for the considered cell, the battery discharge 
rate is a key factor affecting not only the performance variability of the cell, but 
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also the determination of most important random inputs. 
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1. Introduction 

High energy and power densities of Lithium-ion batteries (LIBs) alongside 
their superior safety features have made them the number one energy storage 
device for a wide range of electric devices from cell phones to hybrid-electric 
vehicles and aerospace applications (D l2l El IH. Since the launch of the first 
commercial LIB in 1991 [l5l, significant work has been dedicated to modeling 
[IB 121 mmuni mi, design optimization [|T2l[T3l[T4l . and discovering new materi¬ 
als [[Bl HI 1121 for LIBs. 

Among the proposed models for simulation of LIBs, the most widely used 
one is Newman’s model llT^ [6llT9l, which is based on the porous electrode and 
concentrated solution theories. This model involves material and electrochem¬ 
ical properties, such as porosity of electrodes, solid particle size, and diffusion 
coefficients, which are measured/estimated directly or indirectly via experimental 
techniques [l20l . Measurements of most physical quantities are accompanied by 
uncertainty due to accuracy limitations or natural, cell-to-cell variabilities. For 
example, measurements of the solid particle size in the LIB electrodes typically 
result in distributed values for the particle size [l20l |2T]| rather than a single deter¬ 
ministic value which is mostly used in the LIB simulations. Performing additional 
measurements may lead to a better characterization of the solid particle size, but 
cannot completely eliminate the irreducible uncertainties due to natural variations. 
Hence, a deterministic treatment of the solid particle size may result in predictions 
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which do not agree with the experiments well. Consequently, LIB models whieh 
ineorporate diserete or eontinuous partiele size distributions have been developed 
[|^l^[2^l^[^ . Quantifying the impaet of sueh uneertainties is essential for 
reliable model-based predietions, and is the foeus of the emerging field of Uneer- 
tainty Quantifieation (UQ) in eomputational engineering and seience. In general, 

UQ provides tools for assessing the eredibility of model predietions and faeilitat- 
ing deeision making under uncertainty UQ is also utilized for quanti¬ 

tative validation of simulation models [|29l |30l [3T]| and robust design optimization 
under uneertainty [|^l^[34ll . 

One may represent the uneertain model parameters by random variables/proeesses. 
This is the subjeet of a major elass of UQ approaehes known as probabilistie teeh- 
niques. Among these methods, stoehastie speetral methods UM based on 
polynomial ehaos (PC) expansions [|36l [371 have reeeived speeial attention due to 
their advantages over traditional UQ teehniques sueh as perturbation-based and 
Monte Carlo sampling (MCS) methods. In partieular, under eertain regularity 
eonditions, these sehemes eonverge faster than MCS methods 1(3^ and, unlike 
perturbation methods, are not restrieted to problems with small uneertainty levels 
[l35ll . Stoehastie speetral methods are based on expanding the solution of inter¬ 
est in PC bases. The eoeffieients of these expansions are then eomputed, for in¬ 
stance, intrusively via Galerkin projeetion [l35ll . or non-intrusively via regression 
[l33[40l|4Tl[42l or quadrature integration [|^ . For eomplex systems sueh as LIBs, 
non-intrusive methods are more attraetive than intrusive ones sinee they allow the 
use of simulations as blaek boxes. In other words, there is no need to modify the 
available deterministie solvers when one uses non-intrusive PC expansions. In 
addition, sinee only independent solution realizations are needed, parallel imple- 
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mentation is straightforward. 

Up to date, the majority of the LIB physies-based simulations have treated the 
underlying model parameters deterministieally and ignored the effeets of uneer- 
tainties in the model parameters on the performanee of LIBs. There are a few 
works in the literature whieh have addressed the variability of LIB physies-based 
model parameters. In [l20l . PC expansion is employed within the Newman’s model 
to demonstrate the reduetion in the eell potentials as a result of the uneertainty in 
the partiele size of the anode eleetrode. Effeets of variability in eyeling rate, par- 
tiele size, diffusivity, and eleetrieal eonduetivity of the eathode eleetrode on the 
LIB performanee was examined in [|4^ . where surrogate models are developed 
using the Newman’s model together with teehniques sueh as Kriging, polynomial 
response, and radial-basis neural networks. It was found that the randomness 
in eleetrieal eonduetivity has the minimal effeets on the eell performanee, while 
the impaet of the remaining parameters depend on the the eyeling rate. In fl431l . 
impedanee speetroseopy [|4^ is utilized to eompare the relative importanee of ran¬ 
domness in porosity, partiele size, length and tortuosity of the eathode and sepa¬ 
rator using Nyquist stability eriterion. It was shown that the partiele size has the 
largest impaet on the fluetuations in the impedanee. In WT\ . a redueed order LIB 
model presented in [j9l and a Bayesian inferenee teehnique are used to estimate the 
effeetive kinetie and transport parameters from experimental data. These works 
either employ a redueed order model of the LIB or only eonsider uneertainties in 
a few number of model parameters to avoid high eomputational eosts. 

This work presents a sampling-based PC approaeh to study the effeets of un- 
eertainty in various model parameters on the eell eapaeity, voltage, and eoneen- 
trations of an LIB deseribed by the Newman’s model [|T^ 1^ [T9ll. The proposed 
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PC approach, first introduced in [|40l|, relies on the sparsity of expansion coeffi¬ 
cients to accurately compute the statistics of quantities of interest with a small 
number of battery simulations. Through its application to an LiC6/LiCo02 LIB 
model, we demonstrate that, unlike the previously mentioned works on UQ of 
LIBs, the proposed PC approach is capable of taking into account a large number 
of uncertain parameters. While the proposed PC-based UQ framework is general, 
the results we present are specific to the particular LiC 6 /LiCo 02 cell we consider 
and the the choices of uncertainty models for its parameters. The latter are, as 
much as possible, identified from the reported experiments on identical or similar 
cells/materials. 

Additionally, this UQ framework enables performing a global sensitivity anal¬ 
ysis (SA) to identify the most important uncertain parameters affecting the vari¬ 
ability of the output quantities. Such an analysis may be used toward reducing 
cell-to-cell variations and designing more efficient and targeted quality control 
procedures to reduce the manufacturing cost of LIBs P31l48ll . We also review the 
standard experimental techniques used for measuring the model parameters and 
discuss their impacts on the cell capacity and/or power. 

It is worth highlighting that the present work is focused on modeling and prop¬ 
agation of parametric uncertainties, as opposed to model form uncertainties which 
may consider multiple competing models. 

The remainder of this paper is organized as follows. Section [previews the LIB 
model used in this study. In Section]^ we present our stochastic LIB modeling 
approach which is based on non-intrusive PC expansions. We then continue with 
reviewing a global sensitivity analysis approach via the Sobol’ indices. Standard 
experimental techniques for measuring the model parameters and their underlying 
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distributions are discussed in Section Finally, an LiC6/LiCo02 cell is consid¬ 
ered as our numerical example in Section 


2. LIB governing equations 


An LIB schematic is presented in Fig. During the discharge process, lithium 
ions in the solid particles of the anode diffuse to the particle surface where they 
oxidize into Li+ ions and electrons and transfer to the electrolyte liquid (dein¬ 
tercalation). Electrons flow through the external circuit to the positive electrode. 
Meanwhile, Li’*' ions travel via diffusion and migration through the electrolyte 
and separator to the cathode where they are reduced and diffused into the solid 
particles (intercalation). When the LIB is charged, this process is reversed. 



Macro-scale 


Anode 


Meso-scale 


Separator Cathode 


Micro-scale 


(insertion in cathode 
during discharge) 


(extraction in anode 
during discharge) 


Figure I: Schematic of a full cell LIB. 


The Newman’s LIB model is developed based on porous electrode and con- 
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centrated solution theories HHillll. The governing equations of this model 
for the salt eoneentration in liquid phase c, lithium eoneentration in solid phase 
Cg, liquid phase potential (j)e and solid phase potential are presented in Table 
Symbols used in this table are defined in the nomenelature. The assumptions 
made in this model are summarized as: (i) transport properties are independent 
of temperature and simulations are isothermal; (ii) D, Dg and do not depend 
on the eoneentrations; (iii) volume ehanges within the eell are ignored; (iv) zero 
surfaee eleetrolyte inter-phase (SEI) resistanee is eonsidered; (v) no double layer 
eapaeitive effeets are eonsidered; and (vi) solid partieles are spherieal. 


Table 1: Coupled non-linear governing equations of LIB. 



Governing eqnation 
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(6) 

<1 

n 

II 

o 

: II 

! <1 

n 

II 

II 

o 

Solid phase 

diffusion 

dt r‘^dr\ ^ dr */ 

(7) 

V" Cg 1 j’—0 0 

Vc 1 - 

Liquid 

V(K'=®V<^e) - V(«:g'Vlnc) -f = 0 

(8) 

<1 

d) 

H 

II 

o 

II 

<1 

II 

II 

O 

phase 



03 

II 

II 

O 

potential 




Solid phase 

= 0 

(9) 

0 ^^s|x—L 

potential 



^4^s\x=La ~ ^^s\x=La + Ls ~ 0 


r /0.5Fri\ / 0.5Frj\-\ 

Reaction (^—j - exp (^ - —jj 

kinetics V. = - crT'®(c)“'= 


In the literature, numerieal teehniques sueh as finite differenee Q, finite vol- 
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ume 1491, and finite elements [pQl have been used to solve this system of eoupled 
non-linear equations simultaneously. In order to reduee the eomputational eost of 
the LIB simulations, using the eoneepts of partieular and homogeneous solutions 
to Ordinary Differential Equations (ODEs), Reimers [ISTI suggested a deeoupled 
formulation of Newman’s model, whieh we use in this study. Eor the sake of 


eompleteness, this deeoupling teehnique is presented in Appendix A 


3. UQ via polynomial chaos expansion 

Polynomial ehaos expansion was first introdueed by Wiener in 1938 Il3^ . It 
was reintrodueed to the engineering field in 1991 by Ghanem and Spanos [l35l 
for the problems with Gaussian input uneertainties and later extended to non- 
Gaussian random inputs by using the orthogonal polynomials of the Askey seheme 
(generalized PC expansion) [[53 . PC expansion provides a framework to approx¬ 
imate the solution of a stoehastie system by projeeting it onto a basis of polyno¬ 
mials of the random inputs, whieh we review in the following. 

Eet (fl, T, V) be a eomplete probability spaee, where Vt is the sample set and 
■p is a probability measure on the a—field T. Also assume that the system input 
uneertainty has been diseretized and approximated by random variables, sueh that 
the veetor ^ = (^i, • • • , — )• M'^, d G N, represents the set of independent 

random inputs. We also assume that probability density funetion (PDE) of the 
random variable is denoted by p{ik), while p(^) represents the joint PDE of 
Eet us assume that the finite varianee output quantity of interest (Qol) defined on 
(D, T, V) is denoted by u{^). The truneated PC representation of m(^), denoted 
by m(^), is 

= X] ai'ipiii), (11) 
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where a* are the deterministie eoeffieients and the multivariate PC basis 

funetions. The basis funetions in ( [Tl] ) are generated from 

d 

= Yl'^iki^k), iej^d,p, (12) 

k=l 

where k = 1,... ,(i, are univariate polynomials of degree 4 G No := 

N U {0} orthogonal with respeet to p(4) (see, e.g., Table|5, i.e., 

= j i’iki^k)'ipjt{^k)p{^k)d^k = (13) 

where is the Kroneeker delta and E[-] denotes the mathematieal expeetation 
operator. The multi-index i in ( [IT] ) is i = ( 4 , ■ ■ • , 4 ) ^ d^d,p and the set of 
multi-indiees J4,p is defined by 

d^d,p = {i = ( 4 , • • • ,id) e No : ||i||i 4 p}, (14) 

where || ■ ||i is the h norm and the size of J4,p, henee the number P of PC basis 
funetions of total order not larger than p in dimension d, is given by 


P — \d^d,p\ — 


{p + d)\ 
p\d\ 


(15) 


Due to the orthogonality of the polynomials 4ife('Cfc) and given that the 4 are 
independent, the PC basis funetions are also orthogonal, i.e., E[4i4i] = 


(5ijE[4|]. The truneated PC expansion in (11) eonverges in the mean-square sense 
as p —)• cx) when m(^) has finite varianee and the eoeffieients a* are eomputed 
from the projeetion equation = E[u{-)'ipi{-)]/E[^f\ [[53 . For eonvenienee, 
we also index the PC basis funetions by {1,..., P} so that there is a one-to-one 
eorrespondeneebetween {4j(^)}f=i and 


p,d * 
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Table 2: Correspondence of Wiener-Askey polynomial chaos and probability dis¬ 
tribution of the random variables [15^ . 


p(6) 

Polynomial type 

Support 

Gaussian 

Hermite 

(-CXD,-I-CX)) 

Gamma 

Laguerre 

(0,-1-cx)) 

Beta 

Jacobi 

[a,b] 

Uniform 

Legendre 

[a,b] 


3.1. Non-intrusive polynomial chaos expansion 

The main task in PC-based methods is to compute the coefficients of the so¬ 
lution expansion either intrusively llSfill or non-intrusively [l5^ . In an intrusive 
approach, the governing equations are projected onto the subspace spanned by 
the PC basis via the Galerkin formulation. The final system of equations to be 
solved in an intrusive PC expansion method is P times larger than the size of the 
deterministic counterpart. This approach may require some modifications of the 
existing deterministic solvers, which for complex problems such as LIBs, may 
be difficult and time-consuming to implement. On the other hand, non-intrusive 
methods facilitate the use of existing deterministic solvers and treat them as a 
black box. The first task is to generate a set of N deterministic or random samples 
of denoted by Next, corresponding to these samples, N realizations 

of the output Qol, are computed using an available deterministic 

solver. The last step is solving for the PC coefficients using these realizations. 
Several methods such as least squares regression [[391 , pseudo-spectral colloca¬ 
tion [|29l . Monte Carlo sampling [|38l, and compressive sampling (CS) [l40l have 
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been employed in the literature for this purpose. Onee the PC eoefficients are eom- 
puted, the mean, E[-], and varianee, var[-], of m(^) ean be direetly approximated 
by 


E[u] = ao, 

(16) 

var[fi] = a-. 

(17) 

d,p 


i^O 



In the following, we will review the least squares regression and the eompres- 
sive sampling methods. 

3.1.1. Least squares regression 

The least squares regression teehnique is basieally the regression of the exact 
solution m(^) in the PC bases [l54l . Given the set of samples generated 

randomly, for instance, according to p(^), and the corresponding solution realiza¬ 
tions the discrete representation of ( [IT] ) can be written as 

u, (18) 

where u = G contains the realizations of the Qol, 

is the measurement matrix containing samples of the 
PC basis, and ex = (ai, • • ■ , ap)^ G is the vector of PC coefficients. 

PC coefficients may be approximated by solving the least squares problem 

min II w — ^a||2, (19) 

a 

where || ■ II 2 is the I 2 norm. When is full rank, the solution to ( [T^ is computed 
from the normal equation 

(«''^«')a = (20) 
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In general, a stable solution ck to (20) requires N > P realizations of m(^). In 
[l42l Theorems 2.2 and 3.1], it is shown that for d-dimensional Legendre polyno¬ 


mials of total order p, for the ease of d > p, a stable solution reeovery from (19) 
ean be guaranteed with a number of samples given by 


N > 3PC P\og{P), 


( 21 ) 


where C is an absolute eonstant. This suggests that the number of samples N 
depends linearly on P (up to a logarithmie faetor) for the least squares regression 
method. 

For high dimensional complex problems, such as LIBs, generating N > P 
realizations may be computationally expensive. In such cases, when the solution 
m(^) depends smoothly on the PC coefficients are often sparse, i.e., many of 
them are negligible. In these cases, CS may be employed to compute the coeffi¬ 
cients with N < P realizations [|40l |42l [55l HU |56l |571 |58l ISH . Specifically, the 
importance sampling approach of IHTI ensures an accurate computation of a with 
a number of solution realizations that depends linearly (up to a logarithmic factor 
in P) on the number of dominant coefficients. We next review the CS approach, 
which we use later for the UQ of LIBs. 


3.1.2. Compressive sampling 

Compressive sampling/sensing is an emerging direction in signal processing 
which enables (up to) exact reconstruction of signals admitting sparse represen¬ 
tations with a small number of signal measurements ll60l IMl 1^ . It was first 
introduced to the UQ field in 2011 by Doostan and Owhadi ll4^ where they used 
CS to approximate sparse PC solutions to stochastic PDFs. The main requirement 
in CS methods is to have a sparse solution at the stochastic level, that is a small 
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fraction of PC coefficients in ( [TT] ) are dominant and contribute to the solution 
statisties. The ultimate goal of CS is to approximate the sparse PC eoeffieients ex 
aeeurately and robustly with N < P realizations of m(^). 


With N < P, the underdetermined linear system in (181 is ill-posed and gen¬ 
erally has infinitely many solutions. Sparsity of the PC eoeffieients ex allows a 

I. This ean be aehieved 


regularization of (181 to ensure a well-posed solution 
by solving the Basis Pursuit Denoising (BPDN) problem 

min||Q:||i subject to —it Ih ^ 7- 


( 22 ) 


Minimization of the li norm in (22) promotes sparsity in ex, while the I 2 resid¬ 


ual norm eontrols the aeeuraey of the truneated PC expansion with the toleranee 
7. Several numerieal teehniques are available in the literature to solve the BPDN 
problem Among those, we adopt the Speetral Projeeted Gradient algorithm 
(SPGLl) of [l63l implemented in the SPGLl package for MATLAB 


The accuracy of ex computed from (22) depends on the sample size N and 


the truneation error — 11112 in (22). In general, the truneation error may be 


decreased by increasing the order p of the PC basis, whieh leads to a larger number 
of eoeffieients, P. This, in turn, requires a larger number of samples N to maintain 
the stability of the BPDN problem. The minimum sampling rate depends on the 
type of the PC basis, the sparsity of ex, and the sampling distribution aeeording 
to whieh are generated, and is shown to linearly depend on the number 

of dominant PC eoeffieients iHTl Theorems 3.1 and 4.2]. More preeisely, for d- 
dimensional Legendre polynomials of total order p, where d > p, the number of 


samples N required to guarantee a stable solution recovery from (22) is given by 


N > 3PC S\og{P), 


(23) 
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where S is the number of dominant eoeffieients. In ( [2^ , N depends primarily on 
S and weakly on P through the log(P) term. For high dimensional eomplex prob¬ 
lems, we usually have P S> 5, henee, a eomparison of ( |2T] ) and ( |2^ suggests that 
the CS approaeh requires eonsiderably smaller number of samples in eomparison 
to the least squares regression. 

In praetiee, one may start by approximating a lower order PC expansion when 
N is small and inerease p when a larger number of samples beeome available. 
Another important faetor in the sparse reeonstruetion is the seleetion of the trun- 
eation error toleranee 7. The ideal value for the toleranee is 7 ~ H^'CKexact — 
Sinee the exaet PC eoeffieients c^exact are not known, seleeting larger values 
than ll^'CKexact — wII 2 for 7 deteriorates the aeeuraey of the approximation, while 
smaller ehoiees may result in over-fitting the solution samples and, thus, less aeeu- 
rate results. In the numerieal results of Seetionj^ we employ the eross-validation 
approaeh in @01 Seetion 3.5] to optimally ehoose 7. For the sake of eompleteness, 
this eross-validation approaeh is summarized in Algorithm For more details 
about the CS method, the interested reader is referred to [14011411 . 
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Algorithm 1 Algorithm for cross-validation estimation of 7. 

1 ; Divide the N solution samples to Nr reconstruction and validation sam¬ 
ples. 

2: Choose multiple values for 7^ sueh that the exact truncation error || ^'CKexact — 
1^112 of the reeonstruetion samples is within the range of 7 ^ values. 

3: For eaeh value of 7 ^ solve the BPDN problem ( [22] ) using the Nr reeonstruetion 
samples to eompute ck^- 

4: For eaeh value of 7 ^, eompute the truneation error := — Ur \\2 of 

the Nr validation samples. 

5; Find the minimum value of 7 ^ and and its eorresponding 7 ,, := 7 ^. 



3.2. Global sensitivity analysis 

Identification of the most important random inputs affeeting the variations in 
the cell voltage, capaeity, and concentrations is one of the objectives of this study. 
This is aehieved by performing a global sensitivity analysis (SA) to quantify the 
speeifie effeets of random inputs on the varianee of the Qol. 

Among the available teehniques to perform global SA, we use the Sobol’ in¬ 
dices [[6311 which are widely used due to their generality and accuracy. Sudret 
[[5^ introduced an analytie approaeh to eompute the Sobol’ indiees as a post¬ 
processing of the PC coefficients. Let us assume the PC eoeffieients in ( [TT] ) are 
eomputed. The first order PC-based Sobol’ index Sk, which represents the sole 
effeets of the random input on the variability of m(^), is given by 


Sk='^ af/vaT[u], j 4 = {* e No : 4 > 0,= 0}, (24) 
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where var[M] is given in ini' In eomputing Sk, it is assumed that all random inputs 
exeept are fixed, therefore, Sk does not represent the effeets of the interaetions 
between and other random inputs. In order to quantify the total effeets of the 
random input ^k, ineluding the interaetions between random inputs on the vari¬ 
ability of m(^), one needs to eompute the total PC-based Sobol’ indiees defined 
as 

= {* e : 4 > 0}. (25) 

The smaller S^, the less important random input ^k- For the eases when 
S'J 1, is eonsidered as insignifieant and may be replaeed by its mean value 
without eonsiderable effeets on the variability of m(^). In this study, we employ 
S'[ as a measure to identify the most important random inputs of the LIB model 
eonsidered. 

4. Uncertainty in LIB model parameters 

The model parameters of LIBs are measured experimentally and are aeeom- 
panied by uneertainty due to natural or experimental variability. Some of these 
parameters are measured using eomplex eleetroehemieal teehniques, while others 
are obtained via simple experiments. In the following, we will diseuss a number 
of sueh teehniques. It should be noted that beeause of the limited data available 
in the literature, we eould not avoid making assumptions on the probability distri¬ 
bution for some of the parameters. Additionally, the reported uneertainty models 
are speeifie to the LiC 6 /LiCo 02 cell we consider here. 
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4.1. Porosity, e 

Porosity is defined as the ratio of the pore volume to the bulk volume. Man- 
ufaeturers ehoose the porosity as a trade off between power and energy, i.e., the 
higher the porosity, the higher the power and the lower the eapaeity [|6^ . There 
are several methods to measure the porosity, sueh as the Method of Standard 
Porosimetry (MSP) [|^ . porosity measurement using liquid or gas absorption 
methods aeeording to the Ameriean Soeiety for Testing and Materials (ASTM) 
D-2873 [|M1 . and X-ray tomography [|69l . 

In [143, a uniform distribution [0.28, 032] is eonsidered for the porosity of 
LiCo02 based on experimental data, whieh has ±6.7% (of the mean) variation 
around the mean. DuBeshter et al. [|^ reported ±4.5% variability for poros¬ 
ity of the graphite anode eleetrode, while a higher value of ±12.6% is reported 
for the separator. Henee, in the present study, we assume a uniform distribution 
for the porosity with ±6.7%, ±4.5%, and ±12.6% variation around the mean in 
LiCo02 eathode, LiCg anode, and separator, respeetively. Defining the eoeffieient 
of variation (COV) as the ratio of the standard deviation to the mean, the assumed 
variations translate to COVs of 0.026, 0.073, and 0.038 in anode, separator, and 
eathode, respeetively. 

4.2. Solid particle size, 

The flux of Li+ at the eleetrode-eleetrolyte interfaee is affeeted by the solid 
partiele size in eleetrodes, as it defines the available surfaee area for the reaetion. 
The maximum battery power may be inereased by deereasing the partiele size of 
the eleetrode material and inereasing the surfaee area per volume. 

The partiele size distribution may be measured by a laser diffraetion and seat- 
tering deviee llTOl . X-ray eomputed tomography (XCT) [TTn . or foeused ion beam 
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tomography TF2\ . 

Santhanagopalan and White ll20l considered a normal distribution for a graphite 
anode with the mean and standard deviation of 6.2 /rm and 0.42 /rm, respectively, 
to quantify the effects of random particle size. For the particle size distribution 
of LiCo02, a normal distribution with mean of 7.7 /xm and standard deviation 
of approximately 1.5 /xm is reported in [|2T]| . The corresponding COV values 
are 0.0677 and 0.1948 for the graphite anode and LiCo02 cathode, respectively. 
Nominal values, i.e., mean, of the particle sizes in both electrodes are equal to 
2.0 /xm in our cell. Because of the limited data available in the literature on the 
particle size distribution of our specific electrodes, we use these COVs to find the 
corresponding standard deviations. Based on these COVs, we assume a normal 
distribution with a mean of 2.0 /xm and standard deviation of 0.1354 /xm for the 
graphite anode. For the cathode electrode, the same mean value but a standard 
deviation of 0.3896 /xm is considered. 

Remark 1. When one draws samples from a Gaussian distribution for the solid 
particle size r^, one should assure that all the samples are strictly positive. 
For this purpose, we here employed a truncated Gaussian distribution bounded 
between mean ±3 standard deviation offg. 

4.3. Bruggeman coefficient, brugg 

In porous electrode theory, instead of specifying the exact position and shapes 
of all pores and particles, a volume-averaged formulation is used O. In this 
model, the effective transport properties of the liquid phase, and k®®, are 
obtained using the porosity e and tortuosity r via 
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( 26 ) 


D^« = tR. 

T T 

Tortuosity r is a geometric parameter and depends on the porous electrode struc¬ 
ture [1731 . Although in recent years researchers have attempted to measure r via 
experimental techniques (17311741 [67l [6^ . because of the experimental complex¬ 
ities, T has been commonly used as a model parameter that is calibrated by ex¬ 
periments m- More precisely, r is computed via the well-known Bruggeman 
relation 

T = (27) 

where the Bruggeman exponent brugg is chosen to match numerical results with 
experimental data, and is usually assumed to be 1.5 (1731 . 

Since the Bruggeman relation is widely used in LIB simulations, instead of 
taking r as a random input parameter, we choose the Bruggeman exponent brugg 
to be an uncertain parameter in determining the effective properties. In (IFTII . a uni¬ 
form distribution with ±4.9% of the mean variation around the mean is reported 
for the Bruggeman exponent of the LiCg anode, while ±33.3% is considered for 
the Bruggeman exponent in separator brugg^ [l45ll . For this study, we assume a 
uniform distribution for the Bruggeman exponent with ±5.0% variations around 
the mean in both electrodes and ±20.0% in separator, which corresponds to COVs 
of 0.029 and 0.115, respectively. 

4.4. Li^ transference number, 

When an LIB is discharged, electrolyte salt dissociates into Li+ and PFg ions. 
The portion of the current that is carried by Li+ ions is called the Li+ transference 
number. The higher the Li"*" transport number, the higher the battery power (17511 . 
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Hittorf’s method EH, ac impedance spectroscopy ETII . or pulsed field gradi¬ 
ent NMR (pfg-NMR) technique are used to measure the transference number 
in electrolytes. The experimental error of pfg-NMR technique for binary solu¬ 
tions is estimated to be around 5.0% E9l . Hence, we let the transference number 
uniformly change between 0.345 and 0.381 in this study. 

4.5. Salt diffusion coefficient in the liquid phase, D 

The salt diffusion coefficient D is a measure of the friction forces between 
the ions and the solvents [[80ll . In order to restrict the performance-limiting salt 
concentration gradients, which form in the electrolyte during polarization, it is 
critical to have a high salt diffusion coefficient IfSTlI . High values of D lead to 
higher battery power. Experimental methods such as cyclic voltammetry (CV) 
[[82ll and electrochemical impedance spectroscopy (EIS) ll83ll have been used to 
measure D. 

Salt diffusion coefficient is usually reported as a constant, but in E9ll it is 
assumed to be a function of temperature and salt concentration. In this study, 
we assume a uniform distribution for D with ±10.0% variation around the mean, 
which corresponds to COV of 0.0577. 

4.6. Diffusion coefficient of the solid, Dg 

Diffusion coefficient of the solid phase plays an important role in the perfor¬ 
mance of EIBs since it affects the intercalation flux. Electrochemical techniques 
such as CV, GITT and EIS have been used to determine Dg lIMl . Dg has been 
mostly treated as a constant value in EIB simulations ifT^ . while it has been shown 
that Dg depends on the intercalation level, i.e., the ratio of /cg^rnax ttSH- The 
values of Dg for the same materials, reported by different research groups, may 
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differ by several orders of magnitude. These large differenees may suggest that a 
Fiekian diffusion model whieh is used to ealibrate Dg does not eorreetly deseribe 
the transport of Lithium within the aetive partieles and other transport models 
need to be eonsidered. 

In this study, we assume that Dg does not depend on the interealation level. 
Henee, for Dg, we assume a uniform distribution with the same COV we assumed 
for D. 

4.7. Electronic conductivity of the solid, a 

Capaeity of LIBs may be improved by inereasing the solid phase eleetronie 
eonduetivity by using eonduetive additives in the eleetrode materials If8^ . Two- 
point and four-point probe teehniques have been used to measure the eleetronie 
eonduetivity of eleetrode materials [f87l . Although a depends on temperature and 
state of the eharge [[88l . it is mostly treated as a eonstant in LIB simulations. 

Reported values of a in the literature for LiCe anode are mostly around 100 
S ■ m“^ miElini, while the eleetronie eonduetivity of LiCo 02 eathode is reported 
to be 100 S ■ m“^ 0 and 10 S ■ m“^ In our LIB model, the nominal values of 
a is equal to 100 S ■ m“^ for both eleetrodes. Beeause of the limited available data 
in the literature on the measurement uneertainties in cr, we assume that in both 
eleetrodes it ehanges uniformly with ±10.0% variation. 

4.8. Reaction rate constant, k 

Exehange eurrent density iex is measured at the initial state of the battery ll8^ . 
Using the relation iex = Fk{(ff'^^)^-^{cg^rnax — one ean ealeulate 

the reaetion rate eonstant k for eaeh eleetrode at the initial values of and c. 
Although k shows an Arrhenius type dependenee on temperature and depends on 
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the nature of the eleetrode surfaee [|90ll , it is usually treated as a eonstant in LIB 
simulations. Higher reaetion rate eonstants are favored for Li-ion batteries sinee 
the reaetion is more reversible, and polarization effeets are lower. 

Beeause of the laek of experimental data on the reaetion rate eonstant, we 
assumed that k has a uniform distribution with ±10.0% variation around the mean 
in both eleetrodes. 

A list of random parameters eonsidered in this study is presented in Table 
1^ In addition to the diseussed input parameters, we also took the lengths L of 
eleetrodes and separator to be random in order to study the effeet of uneertainties 
in geometrieal parameters of the battery. It also should be noted that in our battery 
model, the effeetive ionie eonduetivity of the liquid phase is eonsidered to be 
a funetion of liquid eoneentration. Henee, will be automatieally a random 
parameter sinee the porosity is eonsidered to be random. We note that the anodie 
and eathodie transfer eoeffieients (eonsidered to be 0.5 in this study) as well as the 
open eireuit potentials are also subjeeted to uneertainties; however, we treat them 
as deterministie parameters in this study. 
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Table 3: List of random and deterministic LIB inputs used in this study. 


Random Input 

Nominal Value 

Distribution 

rs,a [/um] 

2 

Gaussian, = 2, a = 0.1354 

rs,c [yum] 

2 

Gaussian, /i = 2, a = 0.3896 


0.485 

Uniform, [0.46, 0.51] 


0.724 

Uniform, [0.63, 0.81] 

Cc 

0.385 

Uniform, [0.36, 0.41] 

brugg^ 

4 

Uniform, [3.8, 4.2] 

brugg^ 

4 

Uniform, [3.2, 4.8] 

brugg^ 

4 

Uniform, [3.8, 4.2] 


0.363 

Uniform, [0.345, 0.381] 

D [m^ ■ s“^] 

7.5 X 10-1° 

Uniform, [6.75, 8.25] xl0-i° 

Ds,a [m^ ■ s“^] 

3.9 X 10-1^ 

Uniform, [3.51,4.29] xlO-n 

Ds,c [m^ ■ s“^] 

1 X 10-11 

Uniform, [0.9, 1.1] xlO-n 

aa [S ■ m“i] 

100 

Uniform, [90, 110] 

(Tc [S ■ m“^] 

100 

Uniform, [90, 110] 

ka [m"^ ■ mol ■ s] 

5.03 X 10-11 

Uniform, [4.52, 5.53] xlO-n 

kc [m^ ■ mol ■ s] 

2.334 X 10-11 

Uniform, [ 2.10, 2.56] xlO-n 

La [/im] 

80 

Uniform, [77, 83] 

Ls [/im] 

25 

Uniform, [22, 28] 

Lc [/im] 

88 

Uniform, [85, 91] 
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5. Numerical Example 


In this section, we present a numerical example to demonstrate the applica¬ 
tion of our proposed UQ approach to LIBs. The one-dimensional LiCe/LiCo02 
cell we consider here is studied in [|9l|. We present our results for three different 
discharge rates of 0.25C, 1C, and 4C to study its effects on the propagation of un¬ 
certainty. We note that assuming a constant discharge rate is an idealized scenario. 
Typically, the battery loading varies over the course of the discharge process, de¬ 
pending on an application-dependent usage of the battery. Hence, in addition to 
the previously mentioned uncertain parameters, the discharge rate should also be 
considered as an uncertain input for the LIB model. For the sake of comparison, 
we treat the discharge rate as a deterministic input in this study, while our pro¬ 
posed UQ framework can incorporate random discharge rates with no difficulties, 
perhaps, with a cost of running additional battery simulations. 

For the spatial discretization of the LIB governing equations, we used the 
Finite Difference Method on the non-uniform grids described in [[511 . We also 
performed a mesh convergence analysis to ensure that spatial discretization errors 
are inconsequential. 

Remark 2. Working with a fine mesh and a small time step is advised when one 
employs sampling-based UQ techniques since a coarse discretization may not re¬ 
turn a converged/accurate solution for all samples. In our calculations, we found 
that having 120 non-uniform grid points in each region and 26 non-uniform grid 
points in the spherical solid particles are sufficient to obtain accurate realizations 
for all of the samples. Moreover, constant time steps of 1.0, 0.1, and 0.05 seconds 
were used for 0.25C, 1C, and 4C rates of discharge, respectively. 
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The cell model incorporating the nominal values in Table is referred to as 
the deterministic model, while the stochastic model uses the distributed input pa¬ 
rameters of Table The number of random inputs of our stochastic LIB model 
is d = 19, which may be considered high. We use the CS technique in this study 
to approximate the sparse PC coefficients since, in comparison to least squares 
regression, CS requires smaller number of battery simulations. 

Having at hand a deterministic LIB solver, the next step in sampling-based 
PC expansion is to generate N realizations of model parameters listed in Table 
3. To this end, we assign an independent random variable ^k, k = 1,..., d, to 
each parameter, and consider an appropriate linear transformation of each to 
match the PDF of the corresponding parameter in Table For Gaussian and 
uniform PDFs, we use that are standard Gaussians and uniformly distributed 
between [—1,1], respectively. The independent samples of i.e., are 

used to generate N independent samples of the model parameters, for which the 
LIB model is simulated to obtain N realizations of the output Qols, Me-’niv 
In our numerical experiments, we simulate each battery realization until a cut-off 
potential of 2.8 V is reached. Then using and we solve the 


BPDN problem in (22) to approximate the vector of PC coefficients ck. These in 
turn will be used to compute the statistics of m(^), such as the mean and variance 
given in (|T^ and ([T7|), respectively, as well as the total Sobol’ indices defined in 


(25). 


We find that a PC expansion of order p = 3 and N = 1000 battery simula¬ 
tions are needed to achieve a validation error smaller than 1.0% as specified next. 
Following the cross-validation procedure described in BOl Section 3.5], we divide 
the N = 1000 samples into = 900 reconstruction and = 100 validation 
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samples to estimate the optimum value of 7 in ([22|), using whieh we eompute the 


solution OL from (22). To verify the aeeuraey of the resulting PC expansion, we 
eompute the validation error 


relative error = 


Uvh 


(28) 


where is the veetor of = 1000 additional realizations of Qol (not used in 
eomputing ct) and is the measurement matrix eorresponding to Uy. Stated 


differently, the error in (28) determines the aeeuraey of eonstrueted PC model in 
predieting independent Qol realizations. In praetiee, there is no need to generate 
additional samples for validation and we only used a larger number of validation 
samples to demonstrate the aeeuraey and robustness of our proposed method. 


Remark 3. In practice, we rely on the validation error to estimate the solution 
accuracy and decide to possibly increase the total order of PC expansion p. When 
the validation error does not reduce by increasing the number of numerically gen¬ 
erated LIB samples for a given expansion order p, one may achieve smaller vali¬ 
dations errors by increasing p with a cost of increased number of LIB simulations. 

Algorithm summarizes the main steps in our proposed UQ framework for 
LIBs for a fixed PC expansion order p. 
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Algorithm 2 Summary of the main steps in the proposed UQ framework for LIBs. 
1; Identify uncertain parameters and their probability distributions (Section |^. 


2; while the validation error in (28) is larger than a user-defined threshold do 
3 : Generate N realizations of the uncertain parameters based on their PDFs. 

Perform deterministic LIB simulations to obtain the output Qol (Section 

0 ). 

Evaluate the corresponding realization of the PC basis (Section [3T| ). 
Perform the cross-validation approach in Algorithm[^to optimally estimate 


7 ; Solve the BPDN problem in (22) to approximate the PC coefficients. 

8; Increase N. 

9 ; end while 

10: Compute the statistics of the output Qol via ( [T^ and ( [T7] ) and the total Sobol’ 


indices via (251. 


5.1. Results 

5.1.1. Effects of input uncertainties on cell capacity 

Capacity is defined as the available energy stored in a fully charged LIB and is 
one of the most important factors affecting the battery performance. However, ef¬ 
fects of LIB model uncertainties on capacity estimation are not yet fully explored. 

Lig. demonstrates the effects of LIB model uncertainties in estimating the 
cell capacity. In Lig. |^a), cell capacities as a function of the cell voltage obtained 
from deterministic and stochastic models are presented for 0.25C, 1C and 4C 
rates of galvanostatic discharge. Shaded areas around the mean of stochastic cell 
capacity represent three standard deviation intervals. As it can be seen, noticeable 
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Figure 2: Comparison of stochastic and deterministic battery models for I = 
0.25C, 1C and 4C rates of galvanostatic discharge, (a) Stochastic and determin¬ 
istic cell capacities as a function of the cell voltage. The shaded areas are prob¬ 
ability bounds of three standard deviations around the mean; (b) Relative error 
as a function of the cell voltage; (c) PDFs of the cell capacity; (d) Comparison 
of the relative error obtained by the CS method and the level two SC technique 
on a Clenshaw-Curtis grid [1^ for 1C rate of discharge. Level three SC method 
requires 9976 LIB simulations. 
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deviation from the deterministie eapaeity starts when the eell voltage experienees 
a rapid drop-off. As the diseharge rate is deereased, the onset of this deviation 
happens at higher eell voltages, i.e., at ^ceii 3.7 V for / = 0.25C and (j)ceii ~ 
3.3 V for / = 1C. This plot shows that for low to medium rates of diseharge, a 
deterministie simulation of LIB overestimates the eell eapaeity. For 4C diseharge 
rate, the mean of stoehastie eell eapaeity overlaps with the eell eapaeity predieted 
by the deterministie model. One may suspeet that at high diseharge rates, input 
uneertainties have no signifieant effeet on the eell eapaeity, but the probability 
bounds around the mean of eell voltage for / = 4C suggest otherwise. As it ean 
be seen, the largest standard deviation of the eell eapaeity at the end of diseharge 
is for the 4C diseharge rate. 

Fig. I^b) shows the validation error in approximating the PC eoeffieients 
of the eell eapaeity. With p = 3 and d = 19, the number of PC eoeffieients 
P = 1 ^ 19 , 3 ! = 1540, whieh is larger than N = 1000. The aeeurate solution ap¬ 
proximation in this under-sampled setting is aehieved by minimizing the li norm 
of CK in the BPDN problem ( [22] ). 

A eomparison of deterministie (solid line) and distributed (dashed line) eell 
eapaeities is also presented in Fig. |^e). This plot suggests that the probability 
of aehieving deterministie eell eapaeity when the eell is subjeeted to assumed 
input uneertainties is very small for low to medium diseharge rates, while this 
probability is eonsiderably larger for higher diseharge rates. 

An alternative sampling-based approaeh to eompute the PC eoeffieients is the 
stoehastie eolloeation (SC) method [|92j |94l. The main idea behind the SC 

teehnique is to sample the output Qol at partieular points in the stoehastie spaee 
and then approximate the solution via interpolation or the solution statisties by 
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numerical integration. For high-dimensional problems, sparse tensor products 
first introduced by Smolyak [l95ll . which may be built upon the Clenshaw-Curtis 
abscissas, has been proposed to relax the curse-of-dimensionality associated with 
full tensor products [|9T1 . The accuracy of this method is controlled by the so- 
called level parameter. The number of required LIB simulations increases with 
the magnitude of the level parameter. SC may also be used to compute the PC 
coefficients via the projection equation ai = W[u{-)'rpi{-)]/W[^‘f\. In particular, 
the numerator in this equation is a d-dimensional integral over Vt which may be 
computed using tensor product quadrature or cubature rules. Detailed description 
of these techniques is beyond the scope of this paper and the interested reader is 
referred to the provided references. 

In this study, we employ the SC technique with sparse tensor products based 
on Clenshaw-Curtis abscissas to compare with the proposed CS approach. For the 
level parameter of two, one needs 761 LIB simulations in the SC method. Fig. 
I^d) compares the validation error in approximating the PC coefficients of the 
cell capacity obtained by SC and CS methods for the cell discharged at 1C rate. 
As it can be seen, CS with 761 samples results in the validation errors that are 
smaller (about one order of magnitude) than those reported by the SC technique. 
Increasing the level parameter of the SC method to three, in order to increase the 
accuracy, requires 9976 LIB simulations. A comparison between the level two 
SC and CS approaches confirms the advantage of our proposed UQ framework, 
hence, we did not perform level three SC due to the high computational costs. 
Since we use random samples in CS, unlike in the SC technique, the number of 
samples are not dictated by the method and one may freely choose a minimum 
number of additional LIB simulations to achieve the desired accuracy. 
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Figure 3: Global sensitivity analysis of the cell capacity for various discharge 
rates: (a) / = 0.25C; (b) / = 1C; (c) / = 4C. 
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As discussed in Sectionwe compute the total Sobol’ indices given in (25) 
to identify the most important random inputs. We assume that a random input 
with a maximum total Sobol’ index smaller than 0.01 has no significant effects 
on the variability of the output Qol and may be treated as a deterministic input. 
This selection criteria is problem dependent and may be changed depending on the 
application requirements. Fig. [^shows the total Sobol’ indices associated with the 
cell capacity over a course of discharge for I = 0.25C, 1C and 4C. In this paper, 
we only present significant Sobol’ indices, i.e., those with max(S'^) > 0.01, in the 
global SA plots. Moreover, the legends of the Sobol’ index plots are sorted such 
that the first legend corresponds to the largest S'J, while the last legend represents 
the smallest Sobol’ index. The following observations from Fig. [^are worthwhile 
highlighting: 


• Independent of the discharge rate, porosity e and Bruggeman coefficient 
brugg in anode, separator, and cathode are among the most important ran¬ 
dom inputs. In other words, variation in the the cell capacity is highly af¬ 
fected by uncertainty in tortuosity r. 

• For all three discharge rates. Fig. shows that cxa, cxc, Dg^a, Dg^c, ka, 
kc, and Vg^a are insignificant random inputs and their uncertainties have no 
important effects on the variability of the cell capacity. Consequently, ex¬ 
pensive and accurate quality control measures for these parameters are not 
required when one aims at reducing the variations in cell capacity. 

• As the discharge rate increases, uncertainties in the length of electrodes, 
i.e.. La and Lc, become less important, while the effects of variability in the 
length of separator Lg are more pronounced for / = 4C. 
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• For I = 4C, the diffusion coefficient of the liquid phase D is an impor¬ 
tant random input, while for low to medium rates, its corresponding Sobol’ 
index is smaller than 0.01. 

• By increasing the discharge rate, the number of important random inputs at 
the end of discharge increases; two for / = 0.25C, four for J = 1C, and 
seven for / = 4C. In other words, for high discharge rate LIB applications, 
accurate quality control measures are needed for a larger number of LIB 
parameters. 

• For low discharge rates, variability in the solid particle size has no sig¬ 
nificant effects on the variations of the cell capacity. 

In general. Fig. suggests that in determining the most significant random 
inputs on the variations of the LIB capacity, considering the discharge rate is cru¬ 
cial. Again, we emphasize that these observations may change when more accu¬ 
rate LIB models that account for cell degradations and other physical phenomena 
are employed. 

5.1.2. Statistics in normalized time 

Since the input samples are different, the battery realizations reach the cut-off 
potential of 2.8 V at different times. As an example, for 1C rate of discharge, 
some realizations of the battery reach the cut-off potential before t = 3200 s 
while other realizations need more time, e.g., tmax = 3426 s. This asynchronous 
behavior leads to non-smooth dependence of the cell voltage to random inputs 
and deteriorates the accuracy of PC approximation when the battery approaches 
the end of discharge. 
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(a) (b) 

Figure 4: Effects of input uncertainties on the cell voltage for / = 0.25C, 1C and 
4C rates of discharge, (a) Mean and standard deviation of the cell voltage. The 
shaded areas are probability bounds of three standard deviations around the mean; 
(b) Relative error. 
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To overcome this issue, we introduce an uncertain time scaling t*. In order 
to rescale the deterministic realizations, we assume at t = 0, t* = 100, and we 
set t* = 0 when the battery reaches the cut-off potential of 2.8 V. Hence, at 
t* = 100 the battery is fully charged and at = 0 it is fully discharged. This 
rescaling approach enables us to maintain the accuracy of approximation over the 
entire discharge process without increasing the computational cost or complexity 
of the problem. We note that our definition of t* involves a constraint on the cell 
voltage, such that at = 0 the cell voltage for all realizations is equal to 2.8 V, 
which translates to zero variability for the cell voltage at t* = 0. Our t* definition 
imposes no constraints on other Qol such as solid and liquid phase concentrations. 


5.1.3. Ejfects of input uncertainties on cell voltage and concentrations 

Quantification of the effects of input variations on the cell voltage and concen¬ 
trations over the charge/discharge processes could provide a better understanding 
of the cell behavior. In the following, we present the effects of LIB input uncer¬ 
tainties on the stochastic behavior of the cell voltage fcew, liquid phase concen¬ 
tration c, and solid phase concentration at the surface of the solid particle as 
functions of the normalized time t*. In order to study the variations of c in all 
three main regions of the LIB, we compute c at three different locations in the 
cell; middle of anode, separator and cathode. Similarly, is computed in the 
middle of electrodes. 

Fig. I^a) shows the mean and probability bounds of three standard devia¬ 
tions around the mean of cell voltage as the functions of t* for I = 0.25C, 1C 
and 4C. Although at t* = 100, larger variabilities in fceii correspond to higher dis¬ 
charge rates, fcew experiences its largest standard deviation at 1C rate of discharge. 


35 


total sobol index 
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Figure 5: Global sensitivity analysis of the cell voltage for: (a) / = 0.25C; (b) 
/ = 1C; (c) / = 4C. 


highlighting again the importance of discharge rate on UQ analysis of the LIBs. 
The validation error of the PC solution constructed from N = 1000 samples is 
displayed in Fig. Qb)- In fact, N = 1000 samples were enough to accurately 
approximate c and as well. 

The corresponding total Sobol’ indices of the discharge curves in Fig. |^are 
presented in Fig. from which we highlight that: 

• Similar to the cell capacity, porosity e and Bruggeman coefficient brugg 
in all three regions are among the most important random inputs which 
contribute to the variability of the cell voltage for all three discharge rates. 
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• For all three diseharge rates, unlike the eell capaeity, uneertainties in ka, kc, 
and rs,a eontribute to the eell voltage variations, although, their impaets are 
limited to the onset of diseharge for low to medium diseharge rates. 

• The uneertainty in Ua, cTc, Dg^a, and Dg^c has no signifieant effeets on the 
variability of the eell voltage. 

• As the diseharge rate inereases, random parameters in separator, i.e., e*, 
brugg^, and Lg, eontribute more and more to the eell voltage variations, 
sueh that for 4C rate of diseharge, eg beeomes the most important random 
input affeeting the eell voltage variability. A similar behavior was observed 
in Fig. 

• For I = 4C, variations in the eell voltage are highly affeeted by the un¬ 
eertainties in the diffusion eoeffieient of the liquid phase D. Although, for 
low to medium diseharge rates, D is an insignifieant random input, for high 
diseharge rates, its Sobol’ index S]^ grows eonsiderably. Moreover, unlike 
the eell eapaeity, Li"*’ transferenee number has a high eontribution in the 
eell voltage variability, when the battery is diseharged at 4C rate. 

• Similar to the eell eapaeity, effeets of uneertainties in the length of elee- 
trodes. La and L^, on the variations of 0ceii deereases as the rate of diseharge 
inereases. 

Figs. I^a)-(e) show the mean values of c and their three standard deviation 
bounds at three locations for 0.25C, 1C and 4C rates of discharge. As the LIB 
is discharged at higher rates, larger concentration gradients are developed across 
the cell to balance the migration of anions llT^ . as a result of which the liquid 
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Figure 6: Mean and standard deviation of the liquid phase concentration c in the 
middle of anode, separator and cathode for: (a) I = 0.25C; (b) I = 1C; (c) 4C 
rates of discharge. Mean and standard ^^iation of the solid phase concentration 
^surf middle of anode and cathode for: (d) / = 0.25C; (e) / = 1C; (f) 

4C rates of discharge. The shaded areas are probability bounds of three standard 
deviations around the mean. 
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concentration in the eathode approaehes to zero near the end of diseharge for 4C 
diseharge rate (Fig. |^e)). This may lead to zero eoneentrations in a region at the 
baek of the eathode eleetrode, whieh means the aetive materials may not be uti¬ 
lized further. This limits the eapability of the LIB to be diseharged at higher rates. 
On the other hand, high eoneentrations in the anode electrode may be problematie 
when the lithium salt/solvent system used has a solubility limit. For example, sueh 
a limit is 2100 mol/m^ for Lithium Perehlorate in Propylene Carbonate at room 
temperature [fT^ . Consequently, analysis of the variability of liquid phase eon- 
eentration near the end of diseharge may provide a more aeeurate understanding 
of the rate limiting meehanisms assoeiated with diffusion of lithium ions. As ean 
be observed from this figure, our definition of t* does not impose any eonstraints 
on c; henee, unlike the cell voltage, variability of c is not equal to zero at t* = 0. 

As mentioned in Section it has been shown that a and Dg depend on the 
solid phase eoneentration. Taking into aeeount the variability in the solid phase 
eoneentration due to input variations may help to assess sueh dependeneies more 
aeeurately. A similar analysis is performed for the solid phase eoneentration at 
the surfaee of the solid partiele The mean values of bounded by three 
standard deviations in the middle of eleetrodes for I = 0.25C, 1C and 4C rates of 
diseharge are also presented in Figs. |^d)-(f). It ean be observed that the variability 
of in the middle of the anode does not grow signifieantly as the battery is 
diseharged at low and medium rates, while for the high diseharge rate, standard 
deviation of in both eleetrodes is inereased monotonically over the eourse of 
diseharge. 

Fig. 1^ shows the eorresponding total Sobol’ indiees of c in three regions for 
I = 0.25C, 1C and 4C rates of diseharge. We observe that for low to medium dis- 
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charge rates, the influential random inputs on the variability of c vary eonsiderably 
from one region to another. As an instanee, for / = 0.25C, variations in c in the 
middle of eathode are mainly affeeted by uneertainty in ec, brugg^, and D, while 
in the anode, theses variations are affeeted by the uneertainty in D, Cc, brugg^, e^, 
brugg^, Cs, , Lc, and brugg^. In other words, for low diseharge rates, variability 
of c in the anode is affected by the random inputs assoeiated with the other two 
regions, i.e., separator and eathode, while this is not true for the variation of c in 
the eathode. Additionally, we note that: 

• Similar to the eell eapaeity and voltage, uneertainties in cxa, ac, Ds,a, and 
Dg^c have no signifieant effeets on the variability of c. Moreover, although 
reaetion rate eonstants ka and appear in Figs. |^d), (e), and (h), their 
impaet on the variability of c is small. 

• Despite the diseharge rate in all three regions, uneertainty in D and has 
more pronouneed effeets on the variation of c in both electrodes in eompar- 
ison to separator. 

• 5^ is inereased in the eleetrodes and deereased in the separator as the LIB 
is diseharged at higher rates. 

• Unlike the eell eapaeity and voltage, the effeet of uneertainties in the length 
of the eleetrodes La and on the variations of c inereases as the rate of 
diseharge inereases. 

A global SA of in both eleetrodes is also presented in Fig. We see that 
uneertainty in the solid phase diffusion eoeffleient Dg has no signifieant effeets on 
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the variability of cell capacity, voltage and liquid phase concentration. Surpris¬ 
ingly, Ds may be labeled as an insignificant random input even on the variations 
of in both electrodes for all three discharge rates. Moreover, the uncertainty 
in (Ja and Uc also has no impacts on the variability of independent of the dis¬ 
charge rate. As the LIB experiences higher discharge rates, the total Sobol’ index 
of Vs,a decreases, while for an increase in the total Sobol’ index is observed. 

In summary, we see that in our stochastic LIB example, the uncertainty in a 
and Dg has no significant effects on the variability of the cell capacity, voltage 
and concentrations, hence, may be treated as deterministic inputs. Because of 
this, tight quality control measures are not needed for these parameters. On the 
other hand, e and brugg are the most important random inputs independent of the 
discharge rate, and require tight quality control measures to reduce the LIB cell-to- 
cell variations. Moreover, we observe that the relative contribution of uncertainty 
in the model parameters in the overall performance variability is highly affected 
by the battery discharge rate. 

We emphasize that the results presented in this section apply only to the LiC6/LiCo02 
cell we consider in this study subjected to the input variabilities presented in Table 
1^ Any changes in these uncertainties or cell chemistry/configuration may result in 
different observations. This is also the case if a more accurate LIB model, which 
accounts for additional physical and chemical phenomena, is employed. 

6. Conclusion 

A sampling-based UQ approach was introduced to study the effects of in¬ 
put uncertainties on the performance of the LIBs. The proposed UQ approach is 
based on polynomial chaos expansion framework and hinges on a sparse approx- 
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total sobol index total sobol index total sobol index 



(e) 


(f) 


Figure 8: Sensitivity analysis of the solid phase concentration for / = 0.25C, 1C 
and 4C rates of discharge, (a) in the middle of anode with / = 0.25C; (b) 
^surf middle of cathode with / = 0.25C; (c) in the middle of anode with 
I = 1C; (d) in the middle of cathode with I = 1C; (e) in the middle of 
anode with I = 4C; (f) in the middle of cathode with / = 4C. 
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imation technique to achieve an accurate estimation of solution statistics with a 
small number of LIB simulations. Additionally, the proposed method enables one 
to identify the most important random inputs for any Qol such as capacity, volt¬ 
age, and concentrations by performing a global sensitivity analysis via computing 
the total Sobol’ indices. Such an analysis is helpful in designing more efficient 
and targeted quality control measures, from material selection to cell assembly, 
and reducing the manufacturing cost of the LIBs. 

Performance of the proposed UQ approach was explored through its appli¬ 
cation to an LiC6/LiCo02 LIB discharged at three different rates. It was shown 
that the proposed UQ method can accurately compute the variability in the output 
Qols such as the cell capacity, voltage, and concentrations with a small number of 
battery simulations, 1000 in this example. The global sensitivity analysis results 
corresponding to these Qols showed that the identification of the most important 
uncertain inputs is highly affected by the battery discharge rate. For all three dis¬ 
charge rates, we found that the porosity and Bruggeman coefficient are among the 
most significant uncertain parameters in the performance variability of the exam¬ 
ined LIB. 

We acknowledge that the LIB model we considered in this study suffers from 
many unmodeled phenomena such as side reactions, stresses associated with vol¬ 
ume changes, cell degradation, and temperature dependence, which may affect the 
stochastic behavior of the cell. Since the proposed stochastic LIB model is based 
on sampling, in order to include these phenomena in the model, more accurate LIB 
models can be incorporated without changing the overall UQ framework. This 
may, however, increase the number of required battery simulations and, hence, the 
overall computational cost. Additionally, the accuracy of the uncertainty models 
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we adopted for the LIB input parameters is limited by the available experimental 
data. For some of these inputs, data are highly sparse or non-existent. We, there¬ 
fore, resorted to assumptions in deseribing the uneertainty. When data beeomes 
available, these uneertainty models may be improved aeeordingly. 

The proposed UQ framework in this study was developed for forward prop¬ 
agation of the uncertainties in the LIB simulations. Interesting future research 
directions include using this forward UQ framework along with actual experi¬ 
mental data to infer battery model parameters as well as quantitative validation of 
the model itself. 
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Appendix A. A decoupled formulation for LIBs 


The decoupling procedure starts with combining the potential equations for 
the solid phase Q, liquid phase Q and the Butler-Volmer (BV) equation ( fTO] ) in 
order to first reduce the number of equations to be solved. This is motivated by the 
fact that in each electrode, there are only three tightly coupled quantities, i.e., c, Cg 
and the over-potential rj = (f)s — (l)e — V, which need to be solved simultaneously 
[fSTll . Consequently, the following non-linear equation for the over-potential r] can 
be obtained 


V{a^VT]) = aie 


O.SFnx 

eff 
D 


exp 


0.5F?7\ 
RT ) 


K 


^ + ^vi„c(„) + vi/M 


(A.l) 


subjected to the boundary conditions 

Vr7U=o,L = - \^ + ^Viv) 




x=0,L 


n,D 


^ In c(v) - VV(,,) 


-K 


K 


(A.2) 

(A.3) 


^ — L/a^La-\-L/s 

Here, is the harmonic mean conductivity defined by ^ ^ . So 

far, instead of solving the coupled system of equations in Table one needs to 
solve a system of equations including Eqs. (|^, Q and ( |A.1| ) simultaneously. We 
note that all of the nonlinearities are isolated in Eq. ( A.l| ) for the over-potential. 

The next step of deriving the reformulation is decoupling these three equations 
via the concept of particular and homogeneous solutions to ODEs. Application of 
the implicit finite time differencing to the solid phase diffusion equation (|7]) results 
in the following non-homogeneous ODE with non-homogeneous BCs 
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Cs{r, t + At) - (^Dsr‘^-^Cs{r, t + At) j = Cs{r, t), 


d 




Jvol 


aFD. 


One may write the general solution to Eq. (A.4) as 


(A.4) 

(A.5) 


Cs{r, t + At) = c°{r, t + At) + ‘^ci{r, t + At), (A.6) 

ar 

where c° is the solution to the following non-homogeneous ODE with homoge¬ 
neous BCs 

c:(r, t + At) - (^D,r2^c:(r, t + At)) = c,(r, t), (A.7) 

Vc:|,=o = 0, = 0, (A.8) 

and ci is the solution to the following the homogeneous ODE with non-homogeneous 
BCs 

c^,(r,t +At)-^^(D,r2^ci(r,t +At)) =0, (A.9) 

Vci|,=o = 0, = -^. (A.10) 


A similar approaeh may be employed for deeoupling the solution of the elec¬ 
trolyte phase diffusion equation (|^ from the over-potential solution rj and the 
volumetric pore wall flux j^oi- The only difference here is that j^oi in Eq. is a 
position dependent source term while in Eq. Q it appears on the boundary con¬ 
ditions. Consequently, the general solution to the discretized version of Eq. in 
the time domain can be obtained via 
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c{x, t + At) = c°{x, t + At) + 


1-^H 

F 


jvoi{xo)F{x,xo,t + At)dxo. (A.ll) 


Here in Eq. (A.ll), c° and F are the solutions to the following ODEs both 
subjeeted to zero flux boundary eonditions 


e[c°{x, t + At) — c(x, t)] 

At 


= V[eD"“Vc°(a:,t +At)], 


(A.12) 


eF{x, Xo,t + ^t) 

At 


= V[eZi)®®Vc^(a;, xq, t + At)] + 6{x — xq), 


(A. 13) 


with 6{x — xq) being the Dirae delta funetion. 

At this point, solutions to the solid and eleetrolyte phase diffusion equations 
are decoupled from the over-potential and pore wall flux since in solving Eqs. 
( |A.7[ ), ( |A.9[ ), ( |A.12[ ) and ( |A.13[ ), t] and j^oi are not needed. Moreover, for the 
cases when Dg, and e are constant, analytic solutions are available for Eqs. 


(A.9) and (A. 13) IISTl . c{ and F can also be used to obtain a linearized form for 


the terms including V {rj) and In 0 ( 77 ) on the RHS of Eq. (A.l), respectively. In 


conclusion, one needs to solve the decoupled Eqs. ( |A.7[ ), ( |A.9| ), ( |A.12| ), ( |A.13| ), 
and ( A.l| ) instead of solving the coupled system of non-linear equations in Table 
[^for EIB modeling. In addition, if a constant time step is used, Eqs. ( |A.9 ) and 
( |A.13 1 need to be solved once at the beginning of the simulation. The only major 
approximation in deriving this decoupled formulation is the finite time differenc¬ 
ing which is inevitable in numerical simulations. An extended version of this 
decoupled formulation, which includes SEI resistance and double layer capacitive 
effects coupled with a thermal model, is given in [|96l |97l . Interested reader is 
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referred to ffSTlI for more details on linearizing the over-potential equation dATb, 


grid generation and the time evolution strategy. 
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Nomenclature 

D diffusion coefficient of the liquid 
phase [m“^ ■ s“^] 

Ds diffusion coefficient of the solid 
phase [m“^ ■ s“^] 

F Faraday’s constant = 97484 [C ■ 

mol~^] 

/ total current density across the 

stack [amp ■ m“^] 

L width [m] 

N number of samples 

P number of PC basis functions of 

pth total order in dimension d 

Sk first order Sobol’ index of kth 

random input 

S'J total Sobol’ index of kth random 
input 

T temperature [K] 

V open circuit potential of the ac¬ 

tive material [V] 


sample set 
PC coefficients 

^ vector of input random variables 

5 Kronecker delta or Dirac delta 

function 

e porosity of electrodes and stack 

r] over-potential in electrodes [V] 

7 the truncation error tolerance of 

BPDN problem 

K electronic conductivity of the 

liquid phase [S ■ m“^] 

kd liquid phase diffusional conduc¬ 

tivity [S ■ m“^] 

brugg Bruggeman coefficient 

J^d,p set of pth order multi-indices in 
dimension d 

(pceii cell voltage [V] 

0e Li+ ion potential in liquid phase 

[V] 
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(j)s electron potential in the solid k 

phase [V] 

tjji multivariate PC basis functions p 

p{^) probability measure of random 

variable ^ 

a electronic conductivity of the 

solid phase [S ■ m“^] 

r tortuosity 

a active particle surface area per 

unit volume of electrode [m^ ■ 
m“^] 

c salt concentration in liquid 

phase [mol ■ m“^] 

Cs lithium concentration in solid 

phase [mol ■ m“^] 

lithium concentration in solid 
phase at r = [mol ■ m“^] 

d number of random inputs 

iex exchange current density of an 

electrode reaction [amp ■ m“^] 

jvoi volumetric reaction flux in the 

pore walls [amp ■ m“^] 


reaction rate constant [m"^ ■ mol ■ 
s] 

total order of the PC expansion 

micro-scale distance from the 
center of solid particle [m] 

solid particle size [m] 

t time [s] 

Li+ transference number 

X distance from anode [m] 

Superscripts 

eff effective value 

h harmonic mean 

j unit flux portion 

o internal mixing portion 

Subscripts 

i multi-index 

max maximum 

a anode 

c cathode 
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